# Analysis -------------------------------------------------------------
rm(list=ls())


#install.packages("plm")
library(plm); library(lme4); library(rio)
library(glmmML); 
library(dplyr); library(effects); library(ggplot2)
library(ggeffects)

# Descriptive statistics ---------------------------------------------------------

itanes <- import(file=file.path(here::here(), "data", "processed", "panel_descr.dta"))
itanes$sex <- factor(itanes$sex)
levels(itanes$sex) <- lab.sex
itanes$educ <- factor(itanes$educ)
levels(itanes$educ) <- lab.educ
itanes$polinfo <- factor(itanes$polinfo)
levels(itanes$polinfo) <- lab.polinfo
stargazer(itanes, type = "latex",style = "apsr", digits=3,
          out = file.path(here::here(),"output", "table_summaries.tex"))

# sjmisc::descr(descr_cont, out="browser")
# 
# descr_frq <- descr %>% 
#         dplyr::select(age_cat, female, educ3, vote_ref_mix, lr2, income_fam4cat, area, pid)
# descr_frq$vote_ref_mix <- factor(descr_frq$vote_ref_mix, labels = c("No", "Yes"))
# sjPlot::view_df(descr_frq, show.frq=TRUE, show.prc=TRUE, file="descriptives_freq.html")

# Analysis -------------------------------------------------------------------------

load(file=file.path(here::here(), "data", "processed", "panelWithOmega.Rdata"))
# PTV_PDL: Within-id fixed-effects -----------------------
# Both within-ID and time fixed-effects
ptvpdl_fe_1 <-  plm(ptv_pdl ~ omega.30days.std,
                data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
# Both + controls
ptvpdl_fe_2 <- plm(ptv_pdl ~ omega.30days.std + tg,
               data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
# Both + controls
ptvpdl_fe_3 <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
               data = itanes, model = "within", index = c("id", "wave"), effects="twoways")

stargazer::stargazer(ptvpdl_fe_1, ptvpdl_fe_2, ptvpdl_fe_3, 
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(ptvpdl_fe_1, ptvpdl_fe_2, ptvpdl_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID, PDL", "LR distance with PDL, Std"),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table7_ptvpdl.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))



# Government support (0-10): Within-id fixed-effects -----------------------
# Keep only waves 1 and 2
itanes_gov <- itanes %>% 
        filter(wave %in% c(1,2))

gov_fe_1 <-  plm(gov ~ omega.30days.std,
                data = itanes_gov, model = "within", index = c("id", "wave"), effects="twoways")
gov_fe_2 <- plm(gov ~ omega.30days.std + tg,
               data = itanes_gov, model = "within", index = c("id", "wave"), effects="twoways")
gov_fe_3 <- plm(gov ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
               data = itanes_gov, method = "within", index = c("id", "wave"), effects="twoways")
stargazer::stargazer(gov_fe_1, gov_fe_2, gov_fe_3,
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(gov_fe_1, gov_fe_2, gov_fe_3,
                     covariate.labels = c("Omega 30 days, Std", "PartyID, PDL", "LR distance with PDL, Std"),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table7_gov.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

# BERL: Within-id fixed-effects -----------------------
berl_fe_1 <-  plm(berl ~ omega.30days.std,
                data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
berl_fe_2 <- plm(berl ~ omega.30days.std + tg,
               data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
berl_fe_3 <- plm(berl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
               data = itanes, method = "within", index = c("id", "wave"), effects="twoways")
stargazer::stargazer(berl_fe_1, berl_fe_2, berl_fe_3, 
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(berl_fe_1, berl_fe_2, berl_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID, PDL", "LR distance with PDL, Std"),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table7_berl.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

# LR: Within-id fixed-effects -----------------------
lr_fe_1 <-  plm(lr ~ omega.30days.std,
                data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
lr_fe_2 <- plm(lr ~ omega.30days.std + tg,
               data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
lr_fe_3 <- plm(lr ~ omega.30days.std + tg + pid_pdl.std,
               data = itanes, method = "within", index = c("id", "wave"), effects="twoways")
stargazer::stargazer(lr_fe_1, lr_fe_2, lr_fe_3, 
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(lr_fe_1, lr_fe_2, lr_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", ""),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table7_lr.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))


# All other models:
stargazer::stargazer(gov_fe_1, gov_fe_2, gov_fe_3, berl_fe_1, berl_fe_2, berl_fe_3, lr_fe_1, lr_fe_2, lr_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", "LR distance with PDL, Std "),
                     style = "apsr", type = "html", 
                     out = file.path(here::here(),"output", "table8_AllModels.html"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

stargazer::stargazer(gov_fe_1, gov_fe_2, gov_fe_3, berl_fe_1, berl_fe_2, berl_fe_3, lr_fe_1, lr_fe_2, lr_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", "LR distance with PDL, Std "),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table8_AllModels.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), 
                                    c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                                    c("Tg FE", "Yes", "Yes", "Yes","Yes", "Yes", "Yes","Yes", "Yes", "Yes")))
stargazer::stargazer(gov_fe_3, berl_fe_3, lr_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", "LR distance with PDL, Std "),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table7_AllModels_3only.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "Yes", "Yes", "Yes"), 
                                    c("Year FE", "Yes", "Yes", "Yes"),
                                    c("Tg FE", "Yes", "Yes", "Yes")))

# PTV PD: Within-id fixed-effects -----------------------
# within-ID fixed-effects
fit1.fe <-  plm(ptv_pd ~ omega.30days.std,
                data = itanes, model = "within", index = c("id"), effects="individual")
# Both within-ID and time fixed-effects
fit2.fe <-  plm(ptv_pd ~ omega.30days.std,
                data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
# Both + controls
fit3.fe <- plm(ptv_pd ~ omega.30days.std + pid_pd,
               data = itanes, model = "within", index = c("id", "wave"), effects="twoways")
fit4.fe <- plm(ptv_pd ~ omega.30days.std + pid_pd + tg,
               data = itanes, method = "within", index = c("id", "wave"), effects="twoways")
stargazer::stargazer(fit1.fe, fit2.fe, fit3.fe, fit4.fe, 
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X", "X"), 
                                    c("Year FE", "", "X", "X", "X"),
                                    c("News program FE", "", "", "", "X")))




# Random-effects model (commentare nel testo + Appendix) ---------------------------------------------------------------------------------
ptvpdl_re_1 <- lmer(ptv_pdl ~ omega.30days.std + (1|id),
             data = itanes)
ptvpdl_re_2 <- lmer(ptv_pdl ~ omega.30days.std + tg + (1|id),
             data = itanes)
ptvpdl_re_3 <- lmer(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std + (1|id),
             data = itanes)
stargazer::stargazer(ptvpdl_re_1, ptvpdl_re_2, ptvpdl_re_3, 
                     style = "apsr", type = "text", 
                     omit = c("tg"),
                     add.lines=list(c("Individual RE", "X", "X", "X"), 
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(ptvpdl_re_1, ptvpdl_re_2, ptvpdl_re_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL, Std", "LR distance PDL, Std"),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit = c("tg"),
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table8_ptvpdl_re.tex"),
                     add.lines=list(c("Individual RE", "X", "X", "X"), 
                                    c("Tg FE", "", "X", "X")))


# Suggroup analysis: Frequency of watching tg ---------------------------------------------------------------------------------
itanes_lowfr <- itanes %>% 
        filter(tg_freq==1)
subgroup_tgfreqlow_fe_1 <-  plm(ptv_pdl ~ omega.30days.std,
                data = itanes_lowfr, model = "within", index = c("id", "wave"), effects="twoways")
subgroup_tgfreqlow_fe_2 <-  plm(ptv_pdl ~ omega.30days.std + tg,
                data = itanes_lowfr, model = "within", index = c("id", "wave"), effects="twoways")
subgroup_tgfreqlow_fe_3 <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
               data = itanes_lowfr, model = "within", index = c("id", "wave"), effects="twoways")

stargazer::stargazer(subgroup_tgfreqlow_fe_1, subgroup_tgfreqlow_fe_2, subgroup_tgfreqlow_fe_3, type = "text",
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(subgroup_tgfreqlow_fe_1, subgroup_tgfreqlow_fe_2, subgroup_tgfreqlow_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", ""),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table9_subgr_tgfrqlow.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

itanes_highfr <- itanes %>% 
        filter(tg_freq==3)
subgroup_tgfreqhigh_fe_1 <-  plm(ptv_pdl ~ omega.30days.std,
                                data = itanes_highfr, model = "within", index = c("id", "wave"), effects="twoways")
subgroup_tgfreqhigh_fe_2 <-  plm(ptv_pdl ~ omega.30days.std + tg,
                                data = itanes_highfr, model = "within", index = c("id", "wave"), effects="twoways")
subgroup_tgfreqhigh_fe_3 <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
                               data = itanes_highfr, model = "within", index = c("id", "wave"), effects="twoways")

stargazer::stargazer(subgroup_tgfreqhigh_fe_1, subgroup_tgfreqhigh_fe_2, subgroup_tgfreqhigh_fe_3, 
                     type="text",
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(subgroup_tgfreqhigh_fe_1, subgroup_tgfreqhigh_fe_2, subgroup_tgfreqhigh_fe_3, 
                     covariate.labels = c("Omega 30 days, Std", "PartyID PDL", ""),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table9_subgr_tgfrqhigh.tex"),
                     omit = c("tg"),
                     add.lines=list(c("Individual FE", "X", "X", "X"), 
                                    c("Year FE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

table(itanes$tg_freq)
itanes$tg_freq <- factor(itanes$tg_freq)



# Suggroup analysis: Frequency of watching tg (RE) ---------------------------------------------------------------------------------
itanes$tg_freq <- itanes$tg_freq-8
# -7 "Almost never"
# -6 "once a week"
# -5 "2 week"
# ...3, 4, 5, 6
# 0 "Every day"
subgroup_tgfreq_re_1 <- lmer(ptv_pdl ~ omega.30days.std*tg_freq + (1|id),
             data = itanes)
subgroup_tgfreq_re_2 <- lmer(ptv_pdl ~ omega.30days.std*tg_freq + tg + (1|id),
             data = itanes)
subgroup_tgfreq_re_3 <- lmer(ptv_pdl ~ omega.30days.std*tg_freq + tg + pid_pdl.std + lrdist_pdl.std + (1|id),
             data = itanes)
stargazer::stargazer(subgroup_tgfreq_re_1, subgroup_tgfreq_re_2, subgroup_tgfreq_re_3, 
                     style = "apsr", type = "text", 
                     keep=c("omega.30days.std", "tg_freq", "omega.30days.std:tg_freq", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     add.lines=list(c("Individual RE", "X", "X", "X"), 
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(subgroup_tgfreq_re_1, subgroup_tgfreq_re_2, subgroup_tgfreq_re_3, 
                     covariate.labels = c("Omega (30 days average), Std", "TG frequency","PartyID PDL, Std", "LR distance PDL, Std"),
                     style = "apsr", type = "latex", 
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     keep=c("omega.30days.std", "tg_freq", "omega.30days.std:tg_freq", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table9_subgr_tgfreq_re.tex"),
                     add.lines=list(c("Individual RE", "X", "X", "X"), 
                                    c("Tg FE", "", "X", "X")))

# Plot predicted values
predval <- ggpredict(subgroup_tgfreq_re_3, terms = c("omega.30days.std", "tg_freq"))
ggplot(predval, aes(x=x, y=predicted)) +
        geom_line(size=1, colour="black") + 
        geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha=0.5, fill="grey70") + 
        theme_bw() + scale_color_grey() + 
        facet_wrap(~ group) + 
        ylim(2,5) + 
        xlab("Omega Scores (30-day average)") + ylab("Predicted Probability to vote for PDL")
ggsave(filename = file.path(here::here(), "output", "Figure_X_predprob_tgfreq.tiff"),
       height = 5, width = 7, units = "in", dpi=600,compression="lzw")

# Separate models for tg_freq==3 and tg_freq==1 and 2 (for 0 all would be missing):
hightgfreq <- itanes %>%  filter(tg_freq==3)
lowtgfreq <- itanes %>% filter(tg_freq==1 | tg_freq==2)

fit1.tgfreq <-  plm(ptv_pdl ~ omega.30days.std,
                data = lowtgfreq, 
                model = "within", index = c("id", "wave"), effects="twoways")
fit2.tgfreq <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
               data = lowtgfreq, 
               model = "within", index = c("id", "wave"), effects="twoways")
fit3.tgfreq <-  plm(ptv_pdl ~ omega.30days.std,
                    data = hightgfreq, 
                    model = "within", index = c("id", "wave"), effects="twoways")
fit4.tgfreq <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
                   data = hightgfreq, 
                   model = "within", index = c("id", "wave"), effects="twoways")
stargazer::stargazer(fit1.tgfreq, fit2.tgfreq, fit3.tgfreq, fit4.tgfreq, style = "apsr", type = "text")


# Suggroup analysis: Education ---------------------------------------------------------------------------------
table(itanes$educ)
itanes$educ <- factor(itanes$educ)
levels(itanes$educ) <- c("Primary", "Vocational", "High-School", "Tertiary")

subgroup_educ_re_1 <- lmer(ptv_pdl ~ omega.30days.std*educ + (1|id),
             data = itanes)
subgroup_educ_re_2 <- lmer(ptv_pdl ~ omega.30days.std*educ + tg + (1|id),
             data = itanes)
subgroup_educ_re_3 <- lmer(ptv_pdl ~ omega.30days.std*educ + tg + pid_pdl.std + lrdist_pdl.std + (1|id),
             data = itanes)
stargazer::stargazer(subgroup_educ_re_1, subgroup_educ_re_2, subgroup_educ_re_3,
                     style = "apsr", type = "text",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     keep=c("omega.30days.std", "educ", "omega.30days.std:educ", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     add.lines=list(c("Individual RE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(subgroup_educ_re_1, subgroup_educ_re_2, subgroup_educ_re_3,
                     covariate.labels = c("Omega (30 days average), Std", "PartyID PDL, Std", "LR distance PDL, Std"),
                     style = "apsr", type = "latex",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     keep=c("omega.30days.std", "educ", "omega.30days.std:educ", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table9_subgr_educ_re.tex"),
                     add.lines=list(c("Individual RE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

# Plot predicted values
predval <- ggpredict(subgroup_educ_re_3, terms = c("omega.30days.std", "educ"))
ggplot(predval, aes(x=x, y=predicted)) +
        geom_line(size=1, colour="black") +
        geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha=0.5, fill="grey70") +
        theme_bw() + scale_color_grey() +
        facet_wrap(~ group) +
        ylim(2,5) +
        xlab("Omega Scores (30-day average)") + ylab("Predicted Probability to vote for PDL")
ggsave(filename = file.path(here::here(), "output", "Figure_X_predprob_tgfreq.tiff"),
       height = 5, width = 7, units = "in", dpi=600,compression="lzw")

# Suggroup analysis: Interest in Politics ---------------------------------------------------------------------------------
table(itanes$polint)
itanes$polint <- factor(itanes$polint)
levels(itanes$polint) <- c("Not at all", "Little", "Somewhat", "A lot")

subgroup_polint_re_1 <- lmer(ptv_pdl ~ omega.30days.std*polint + (1|id),
                           data = itanes)
subgroup_polint_re_2 <- lmer(ptv_pdl ~ omega.30days.std*polint + tg + (1|id),
                           data = itanes)
subgroup_polint_re_3 <- lmer(ptv_pdl ~ omega.30days.std*polint + tg + pid_pdl.std + lrdist_pdl.std + (1|id),
                           data = itanes)
stargazer::stargazer(subgroup_polint_re_1, subgroup_polint_re_2, subgroup_polint_re_3,
                     style = "apsr", type = "text",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     keep=c("omega.30days.std", "polint", "omega.30days.std:polint", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     add.lines=list(c("Individual RE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))
stargazer::stargazer(subgroup_polint_re_1, subgroup_polint_re_2, subgroup_polint_re_3,
                     covariate.labels = c("Omega (30 days average), Std", "PartyID PDL, Std", "LR distance PDL, Std"),
                     style = "apsr", type = "latex",
                     star.cutoffs = c(0.05, 0.01, 0.001),
                     digits = 3,
                     keep=c("omega.30days.std", "polint", "omega.30days.std:polint", "pid_pdl.std", "lrdist_pdl.std", "Constant"),
                     omit.stat = c("adj.rsq"),
                     out = file.path(here::here(),"output", "table9_subgr_polint_re.tex"),
                     add.lines=list(c("Individual RE", "X", "X", "X"),
                                    c("Tg FE", "", "X", "X")))

# Plot predicted values
predval <- ggpredict(subgroup_polint_re_3, terms = c("omega.30days.std", "polint"))
ggplot(predval, aes(x=x, y=predicted)) +
        geom_line(size=1, colour="black") +
        geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha=0.5, fill="grey70") +
        theme_bw() + scale_color_grey() +
        facet_wrap(~ group) +
        ylim(2,5) +
        xlab("Omega Scores (30-day average)") + ylab("Predicted Probability to vote for PDL")
ggsave(filename = file.path(here::here(), "output", "Figure_X_predprob_tgfreq.tiff"),
       height = 5, width = 7, units = "in", dpi=600,compression="lzw")





# # Separate models for educ==3 and educ==0, 1 and 2 :
# higheduc <- itanes %>%  filter(educ==3)
# loweduc <- itanes %>% filter(educ==0| educ==1 | educ==2)
# 
# fit1.educ <-  plm(ptv_pdl ~ omega.30days.std,
#                     data = loweduc, 
#                     model = "within", index = c("id", "wave"), effects="twoways")
# fit2.educ <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
#                    data = loweduc, 
#                    model = "within", index = c("id", "wave"), effects="twoways")
# fit3.educ <-  plm(ptv_pdl ~ omega.30days.std,
#                     data = higheduc, 
#                     model = "within", index = c("id", "wave"), effects="twoways")
# fit4.educ <- plm(ptv_pdl ~ omega.30days.std + tg + pid_pdl.std + lrdist_pdl.std,
#                    data = higheduc, 
#                    model = "within", index = c("id", "wave"), effects="twoways")
# stargazer::stargazer(fit1.educ, fit2.educ, fit3.educ, fit4.educ, style = "apsr", type = "text")